# Set the working directory and tool paths on your local computer.
WD <- '/Users/Alec/motrpac/20210826_pass1c-umich'
# Set the gsutil path
knitr::opts_knit$set(root.dir=WD)
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(cache = FALSE)write_css = FALSE
if(write_css){
writeLines("td, th { padding : 3px } th { background-color:white; color:black; border:1px solid black; text-align:center } td {color:black; border:1px solid black; word-wrap:break-word; white-space:nowrap; overflow: hidden; text-overflow: ellipsis; max-width:300px; text-align:left}", con=file.path(normalizePath(WD), "style/style.css"))
}################################################################################
##### Resources and Dependencies ###############################################
################################################################################
# Whether to knit document and display data
knit_time = TRUE
# Load dependencies
pacs...man <- c("tidyverse","kableExtra","devtools","MotrpacBicQC","impute","glue",
"rethinking")
for(pac in pacs...man){
suppressWarnings(suppressPackageStartupMessages(library(pac, character.only = TRUE)))
}
#browseVignettes("MotrpacBicQC")
############################################################
##### Functions ############################################
############################################################
# Name functions
select <- dplyr::select
map <- purrr::map
desc <- dplyr::desc
arrange <- dplyr::arrange
melt <- reshape2::melt
mutate <- dplyr::mutate
glue <- glue::glue
# Global options
options(dplyr.print_max = 100)
options(scipen=10000)
# Colors
redblue100<-rgb(read.table(paste0(WD,'/colors/redblue100.txt'),sep='\t',row.names=1,header=T))ds <- 'pass1a'
site <- 'umich'
tech <- 'rppos'
TIS <- 'PLA'
tis <- 'pla'# Phenotype Data (1A)
#########################
pheno_file <- glue("{WD}/data/20201021_pass1a-06-pheno-viallabel_steep.txt")
pheno_df1a <- read.csv(pheno_file, header = T, sep = '\t')
pheno_df1a <- pheno_df1a %>%
mutate(tissue = case_when(Specimen.Processing.sampletypedescription == 'Brown Adipose' ~ 'badi',
Specimen.Processing.sampletypedescription == 'EDTA Plasma' ~ 'pla',
Specimen.Processing.sampletypedescription == 'Gastrocnemius' ~ 'gas',
Specimen.Processing.sampletypedescription == 'Heart' ~ 'hrt',
Specimen.Processing.sampletypedescription == 'Kidney' ~ 'kid',
Specimen.Processing.sampletypedescription == 'Liver' ~ 'liv',
Specimen.Processing.sampletypedescription == 'Lung' ~ 'lun',
Specimen.Processing.sampletypedescription == 'White Adipose' ~ 'wadi',
Specimen.Processing.sampletypedescription == 'PaxGene RNA' ~ 'pax',
Specimen.Processing.sampletypedescription == 'Hippocampus' ~ 'hip',
Specimen.Processing.sampletypedescription == 'Cortex' ~ 'cor',
Specimen.Processing.sampletypedescription == 'Hypothalamus' ~ 'hyp',
Specimen.Processing.sampletypedescription == 'Vastus Lateralis' ~ 'vas',
Specimen.Processing.sampletypedescription == 'Tibia' ~ 'tib',
Specimen.Processing.sampletypedescription == 'Aorta' ~ 'aor',
Specimen.Processing.sampletypedescription == 'Small Intestine' ~ 'sma',
Specimen.Processing.sampletypedescription == 'Adrenals' ~ 'adr',
Specimen.Processing.sampletypedescription == 'Colon' ~ 'col',
Specimen.Processing.sampletypedescription == 'Spleen' ~ 'spl',
Specimen.Processing.sampletypedescription == 'Testes' ~ 'tes',
Specimen.Processing.sampletypedescription == 'Ovaries' ~ 'ova'))
pheno_df1a$viallabel <- as.character(pheno_df1a$viallabel)Phenotypic Data Available: Kidney, Heart, Lung, Gastroc, Liver, Plasma, White Adipose Abundance Data not available for white adipose
# #created in 20210910_pass1a-umich-sample-annotation_steep.Rmd
# order_file <- glue("{WD}/data/20200910_pass1a1c-sample-order_steep.txt")
# sample.order<-read.delim(order_file,header=T, sep="\t")
# Load the prior pass1a data (takes a few minutes)
pass1a_nested_file <- glue("{WD}/../20200915_metabolomics-pass1a/data/20201010_pass1a-metabolomics-countdata-nested_steep.rds")
pass1a_df <- readRDS(pass1a_nested_file)
# Sample Order
pla_rppos_pass1a.0.order <- pass1a_df %>%
filter(STUDY_INSTITUTE == 'University of Michigan') %>%
filter(NAMED == 'named') %>%
filter(DATASET == 'untargeted') %>%
filter(METAB_FAMILY == 'rppos') %>%
filter(TISSUE == 'plasma') %>%
select(SAMPLE_DATA) %>%
unnest(SAMPLE_DATA) %>%
left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)
rm(pass1a_df)# UM-rppos
load(file = glue("{WD}/data/UM-rppos/UM_rppos.0.RData"))
pla1a.0 <- pla_rppos_pass1a.0# Remove internal standards
is <- colnames(pla1a.0)[grepl("istd", colnames(pla1a.0), ignore.case = TRUE)]
# Subset matrix
pla1a.0 <- pla1a.0[, !colnames(pla1a.0) %in% is]# Collect the distances from reference samples
tmp.iter1a <- pla_rppos_pass1a.0.order %>%
mutate(reference = ifelse(str_sub(sample_id,1,1) == '9', 0, 1)) %>%
mutate(drift = ifelse(sample_type == 'QC-DriftCorrection', 1, 0)) %>%
arrange(sample_order)
tmp.iter1a$right_p <- 0
for(i in 1:nrow(tmp.iter1a)){
# set t
if(tmp.iter1a[i,'reference'] == 1){
t = 0
}else if(tmp.iter1a[i,'reference'] == 0){
t = t + 1
}
tmp.iter1a[i,"right_p"] <- t
}
t=0
tmp.iter1a$right_p_d <- 0
for(i in 1:nrow(tmp.iter1a)){
# set t
if(tmp.iter1a[i,'drift'] == 1){
t = 0
}else if(tmp.iter1a[i,'drift'] == 0){
t = t + 1
}
tmp.iter1a[i,"right_p_d"] <- t
}
t=0
tmp.iter1a <- tmp.iter1a %>%
arrange(desc(sample_order))
tmp.iter1a$left_p <- 0
for(i in 1:nrow(tmp.iter1a)){
# set t
if(tmp.iter1a[i,'reference'] == 1){
t = 0
}else if(tmp.iter1a[i,'reference'] == 0){
t = t + 1
}
tmp.iter1a[i,"left_p"] <- t
}
t=0
tmp.iter1a$left_p_d <- 0
for(i in 1:nrow(tmp.iter1a)){
# set t
if(tmp.iter1a[i,'drift'] == 1){
t = 0
}else if(tmp.iter1a[i,'drift'] == 0){
t = t + 1
}
tmp.iter1a[i,"left_p_d"] <- t
}
t=0
tmp.iter1a <- tmp.iter1a %>%
rowwise() %>%
mutate(min_p = min(c(left_p,right_p) ,na.rm = TRUE)) %>%
mutate(sum_p = sum(c(left_p,right_p) ,na.rm = TRUE)) %>%
mutate(min_p_d = min(c(left_p_d,right_p_d) ,na.rm = TRUE)) %>%
mutate(sum_p_d = sum(c(left_p_d,right_p_d) ,na.rm = TRUE)) %>%
arrange(sample_order)
tmp.join1a <- tmp.iter1a %>%
select(sample_id, sample_order, left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d)
# Collect the sample order for test+ref samples
tmp.ref1a <- pla_rppos_pass1a.0.order %>%
arrange(sample_order) %>%
mutate(control = ifelse(grepl('Control', Key.anirandgroup), 1, 0)) %>%
mutate(drift = ifelse(grepl('Drift', sample_type), 1, 0)) %>%
mutate(reference = ifelse(str_sub(sample_id,1,1) == '9', 0, 1)) %>%
mutate(time = case_when(grepl('IPE', Key.anirandgroup) ~ 0,
grepl('0 hr', Key.anirandgroup) ~ 0,
grepl('0.5 hr', Key.anirandgroup) ~ 0.5,
grepl('1 hr', Key.anirandgroup) ~ 1,
grepl('4 hr', Key.anirandgroup) ~ 4,
grepl('7 hr', Key.anirandgroup) ~ 7,
grepl('24 hr', Key.anirandgroup) ~ 24,
grepl('48 hr', Key.anirandgroup) ~ 48)) %>%
left_join(y = tmp.join1a) %>%
select(sample_id,sample_order, Registration.sex, control, time, drift, reference,
left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d)
N2 <- tmp.ref1a[,1] %>% unlist()
miss1 <- N2[!N2 %in% row.names(pla1a.0)] # Verify all samples are in the pla1a.0 file
miss1## named character(0)
# sample.order %>% filter(sample_id == "90750016606")
print('Vice Versa:')## [1] "Vice Versa:"
miss2 <- row.names(pla1a.0)[!row.names(pla1a.0) %in% N2]
miss2## character(0)
N2 <- N2[!N2 %in% c(miss1,miss2)] # TODO: investigate why samples are missing, continue for now
# Reorder pla1a.0 by run order
pla1a.0 <- pla1a.0[N2,]
all(N2 == row.names(pla1a.0)) # Must be true## [1] TRUE
tmp.ref1a <- tmp.ref1a %>%
filter(sample_id %in% N2) %>%
select(sample_order, Registration.sex, control, time, drift, reference,
left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d) %>% as.matrix()
# Collect the sample order for test samples
options(digits = 14)
tmp.sample1a <- pla_rppos_pass1a.0.order %>%
filter(substr(sample_id, 1, 1) == '9') %>%
arrange(sample_order) %>%
mutate(control = ifelse(grepl('Control', Key.anirandgroup), 1, 0)) %>%
mutate(time = case_when(grepl('IPE', Key.anirandgroup) ~ 0,
grepl('0 hr', Key.anirandgroup) ~ 0,
grepl('0.5 hr', Key.anirandgroup) ~ 0.5,
grepl('1 hr', Key.anirandgroup) ~ 1,
grepl('4 hr', Key.anirandgroup) ~ 4,
grepl('7 hr', Key.anirandgroup) ~ 7,
grepl('24 hr', Key.anirandgroup) ~ 24,
grepl('48 hr', Key.anirandgroup) ~ 48)) %>%
left_join(y = tmp.join1a) %>%
mutate(sample_id = str_replace_all(sample_id, pattern = '-', replacement = '')) %>%
mutate(sample_id = as.numeric(sample_id)) %>%
select(sample_id, sample_order, Registration.sex, control, time,
left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d) %>%
as.matrix()
N1 <- row.names(pla1a.0)[substr(row.names(pla1a.0), 1, 1) == '9']
tmp.sample1a <- tmp.sample1a[tmp.sample1a[,1] %in% N1,]
# Reorder pla1a.0.nr by run order
pla1a.0.nr <- pla1a.0[N1,]
tmp.sample1a[,1][!tmp.sample1a[,1] %in% row.names(pla1a.0.nr)] # Verify all samples are in the pla1a.0 file## numeric(0)
all(as.character(tmp.sample1a[,1]) == row.names(pla1a.0.nr))## [1] TRUE
# If out of order, command below will ensure pla1a.0.nr in run order
#pla1a.0.nr <- pla1a.0.nr[as.character(tmp.sample1a[,1]),]NR <- dim(pla1a.0)[1]
P <- dim(pla1a.0)[2]
dim(pla1a.0)## [1] 101 170
N <- dim(pla1a.0.nr)[1]
dim(pla1a.0.nr)## [1] 77 170
confirmed: no negative or zero values
min(pla1a.0,na.rm=T)## [1] 59.47415952
Blank reference samples at the beginning and end
pla1a.0.f.c0<-apply(pla1a.0,1,function(x) sum(is.na(x)))
plot(pla1a.0.f.c0, ylim = c(0,P))No blank test samples
pla1a.0.nr.f.c0<-apply(pla1a.0.nr,1,function(x) sum(is.na(x)))
plot(pla1a.0.nr.f.c0, ylim = c(0,P))pla1a.0.f.c0<-apply(pla1a.0,2,function(x) sum(is.na(x)))
plot(pla1a.0.f.c0, ylim = c(0,NR))pla1a.0.nr.f.c0<-apply(pla1a.0.nr,2,function(x) sum(is.na(x)))
plot(pla1a.0.nr.f.c0, ylim = c(0,N)); abline(h = 10)Remove high missing features
rm_n <- sum(pla1a.0.nr.f.c0>=10)
pla1a.0 <- pla1a.0[,pla1a.0.nr.f.c0<10]
pla1a.0.nr <- pla1a.0.nr[,pla1a.0.nr.f.c0<10]
dim(pla1a.0)## [1] 101 170
dim(pla1a.0.nr)## [1] 77 170
pla1a.0.1 <- log(pla1a.0, 2)
pla1a.0.nr1 <- log(pla1a.0.nr, 2)sum(is.na(pla1a.0.1))## [1] 0
sum(is.na(pla1a.0.nr1))## [1] 0
feature_impute <- apply(is.na(pla1a.0.nr1),2,sum)[apply(is.na(pla1a.0.nr1),2,sum) > 0]pla1a.0.nr1.f.c0<-apply(pla1a.0.nr1,2,function(x) sum(is.na(x)))
plot(pla1a.0.nr1.f.c0, ylim = c(0,NR))# NR
# N
# P
neg_vals <- 0
zero_vals <- 0
feature_na_filter <- 10
P1 <- dim(pla1a.0.nr)[2]
NR1 <- dim(pla1a.0)[1]
N1 <- dim(pla1a.0.nr)[1]
na_vals_impute <- sum(is.na(pla1a.0.nr1))
knn_k <- 10
feature_impute## named integer(0)
if(na_vals_impute > 0){
glue("Features & Values to impute:")
feature_impute
pla1a.0.nr2<-impute.knn(pla1a.0.nr1,k=10)$data
#view the features before and after imputation
par(mfrow=c(2,1),bg="black")
image(as.matrix(pla1a.0.nr1[,pla1a.0.nr1.f.c0>0]),col=redblue100,axes=F)
image(as.matrix(pla1a.0.nr2[, pla1a.0.nr1.f.c0>0]),col=redblue100,axes=F)
par(mfrow=c(1,1) ,bg="white")
glue("Verified no missing values: {sum(is.na(pla1a.0.nr2))}") #verified 0
}else{
print('No missing values to impute')
pla1a.0.nr2 <- pla1a.0.nr1
}## [1] "No missing values to impute"
# feature_impute2 <- apply(is.na(pla1a.0.1),2,sum)[apply(is.na(pla1a.0.1),2,sum) > 0]
# glue("Features & Values to impute:")
# feature_impute2
# pla1a.0.2<-impute.knn(pla1a.0.1,k=10)$data
# #view the features before and after imputation
# par(mfrow=c(2,1),bg="black")
# image(as.matrix(pla1a.0.1[,pla1a.0.nr1.f.c0>0]),col=redblue100,axes=F)
# image(as.matrix(pla1a.0.2[,pla1a.0.nr1.f.c0>0]),col=redblue100,axes=F)
# par(mfrow=c(1,1) ,bg="white")
# glue("Verified no missing values: {sum(is.na(pla1a.0.2))}") #verified 0a <- 0.92
b <- 1
x <- tmp.ref1a[,1]
sidebar <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
sidebar <- cbind(sidebar,sidebar,sidebar)
cor.tmp<-cor(t(pla1a.0.1),method="spearman",use="pairwise.complete.obs") #includes ref samples
dim(cor.tmp)## [1] 101 101
image(cbind(cor.tmp,sidebar),
col=redblue100,axes=FALSE,zlim=c(0.92,1), main=glue("{TIS}2-{NR1}, run-order, z=(0.92,1)"), asp = 1)if(knit_time){
pla_rppos_pass1a.0.order %>%
arrange(sample_order) %>%
select(sample_id, sample_type, sample_order) %>%
knitr::kable(format = "html") %>%
scroll_box(width = "100%", height = "400px")
}| sample_id | sample_type | sample_order |
|---|---|---|
| CS00000MP-1 | QC-DriftCorrection | 3 |
| R00CHRPL2-8 | QC-InternalStandard | 4 |
| CSMR80015-1 | QC-Reference | 5 |
| CSMR80014-1 | QC-Reference | 6 |
| CS00000MP-2 | QC-DriftCorrection | 7 |
| CS000QCMP-1 | QC-Pooled | 8 |
| 90109013104 | Sample | 9 |
| 90121013104 | Sample | 10 |
| 90156013104 | Sample | 11 |
| 90031013104 | Sample | 12 |
| 90127013104 | Sample | 13 |
| 90133013104 | Sample | 14 |
| 90044013104 | Sample | 15 |
| 90025013104 | Sample | 16 |
| 90138013104 | Sample | 17 |
| CS00000MP-3 | QC-DriftCorrection | 18 |
| 90159013104 | Sample | 19 |
| 90148013104 | Sample | 20 |
| 90114013104 | Sample | 21 |
| 90041013104 | Sample | 22 |
| 90005013104 | Sample | 23 |
| 90017013104 | Sample | 24 |
| 90050013104 | Sample | 25 |
| 90037013104 | Sample | 26 |
| 90142013104 | Sample | 27 |
| CS00000MP-4 | QC-DriftCorrection | 28 |
| CS000QCMP-2 | QC-Pooled | 29 |
| 90010013104 | Sample | 30 |
| 90027013104 | Sample | 31 |
| 90009013104 | Sample | 32 |
| 90123013104 | Sample | 33 |
| 90034013104 | Sample | 34 |
| 90124013104 | Sample | 35 |
| 90008013104 | Sample | 37 |
| 90117013104 | Sample | 38 |
| CS00000MP-5 | QC-DriftCorrection | 39 |
| 90007013104 | Sample | 40 |
| 90112013104 | Sample | 41 |
| 90028013104 | Sample | 42 |
| 90011013104 | Sample | 43 |
| 90026013104 | Sample | 44 |
| 90141013104 | Sample | 45 |
| 90126013104 | Sample | 46 |
| 90115013104 | Sample | 47 |
| 90134013104 | Sample | 48 |
| CS00000MP-6 | QC-DriftCorrection | 49 |
| CS000QCMP-3 | QC-Pooled | 50 |
| 90042013104 | Sample | 51 |
| 90119013104 | Sample | 52 |
| 90023013104 | Sample | 53 |
| 90129013104 | Sample | 54 |
| 90029013104 | Sample | 55 |
| 90020013104 | Sample | 56 |
| 90046013104 | Sample | 57 |
| 90122013104 | Sample | 58 |
| 90039013104 | Sample | 59 |
| CS00000MP-7 | QC-DriftCorrection | 60 |
| 90157013104 | Sample | 61 |
| 90048013104 | Sample | 62 |
| 90032013104 | Sample | 63 |
| 90018013104 | Sample | 64 |
| 90160013104 | Sample | 65 |
| 90155013104 | Sample | 66 |
| 90145013104 | Sample | 67 |
| 90140013104 | Sample | 68 |
| 90033013104 | Sample | 69 |
| CS00000MP-8 | QC-DriftCorrection | 70 |
| CS000QCMP-4 | QC-Pooled | 71 |
| 90110013104 | Sample | 72 |
| 90143013104 | Sample | 73 |
| 90147013104 | Sample | 74 |
| 90038013104 | Sample | 75 |
| 90139013104 | Sample | 76 |
| 90150013104 | Sample | 77 |
| 90053013104 | Sample | 78 |
| 90014013104 | Sample | 79 |
| 90052013104 | Sample | 80 |
| CS00000MP-9 | QC-DriftCorrection | 81 |
| 90146013104 | Sample | 82 |
| 90040013104 | Sample | 83 |
| 90118013104 | Sample | 84 |
| 90120013104 | Sample | 85 |
| 90128013104 | Sample | 86 |
| 90043013104 | Sample | 87 |
| 90136013104 | Sample | 88 |
| 90013013104 | Sample | 89 |
| 90135013104 | Sample | 90 |
| CS00000MP-10 | QC-DriftCorrection | 91 |
| CS000QCMP-5 | QC-Pooled | 92 |
| 90015013104 | Sample | 93 |
| 90047013104 | Sample | 94 |
| 90152013104 | Sample | 95 |
| 90144013104 | Sample | 96 |
| 90012013104 | Sample | 97 |
| 90045013104 | Sample | 98 |
| CS000QCMP-6 | QC-Pooled | 99 |
| CS00000MP-11 | QC-DriftCorrection | 100 |
| CSMR80015-2 | QC-Reference | 101 |
| CSMR80014-2 | QC-Reference | 102 |
| R00CHRPL1-3 | QC-InternalStandard | 103 |
| CS00000MP-12 | QC-DriftCorrection | 106 |
a <- 0.92
b <- 1
x <- tmp.sample1a[,2]
sidebar <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
sidebar <- cbind(sidebar,sidebar,sidebar)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
glue("Verified {N} test samples: {dim(cor.tmp)}") #verified, =78## Verified 77 test samples: 77
## Verified 77 test samples: 77
image(
cbind(cor.tmp[order(tmp.sample1a[,2]),],sidebar),
col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Run-Order, z=(0.92,1)"), asp = 1)plot(apply(cor.tmp[order(tmp.sample1a[,2]),order(tmp.sample1a[,2])],1,median))cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
dim(cor.tmp) #verified, =78## [1] 77 77
plot(apply(cor.tmp,1,median)); abline(h=0.96, col = 'blue')# Determine which samples are outliers
o.n <- c(1:N1)[apply(cor.tmp,1,median)<0.96]
o.s <- colnames(cor.tmp)[apply(cor.tmp,1,median)<0.96]
glue("Outlier Samples: {length(o.s)}")## Outlier Samples: 1
o.s## [1] "90122013104"
if(knit_time){
pla_rppos_pass1a.0.order %>%
filter(sample_id %in% o.s) %>%
knitr::kable(format = "html") %>%
scroll_box(width = "100%", height = "100%")
}| sample_id | sample_type | sample_order | Key.anirandgroup | Registration.sex |
|---|---|---|---|---|
| 90122013104 | Sample | 58 | Exercise - 4 hr | 1 |
o.f <- 0.96The one sample that could be classified as an outlier seems to have a similar expression signature of is male/female contemporaries
# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,3]
sex.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
x <- tmp.sample1a[,4]
control.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
x <- as.factor(tmp.sample1a[,5]) %>% as.numeric()
time.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
sidebar<-cbind(sex.type,sex.type,sex.type,
control.type,control.type,control.type,
time.type,time.type,time.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
cbind(
cor.tmp[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]),
order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5])],
sidebar[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]),]),
col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N}, Sex-Control-Time, z=(0.92,1)"), asp = 1)plot(apply(cor.tmp[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]),
order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5])],
1,median))# Add color sidebar
a <- 0.92
b <- 1
x <- is.na(tmp.ref1a[,2]) %>% as.integer()
sample.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
sidebar<-cbind(sample.type,sample.type,sample.type,sample.type,sample.type)
cor.tmp<-cor(t(pla1a.0.1),method="spearman",use="pairwise.complete.obs") #includes ref samples
image(
cbind(
cor.tmp, sidebar),
col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2, Run-Order (+Ref-type), z=(0.92,1)"), asp = 1)run_var <- 0
outlier_sample_n <- length(o.s)
outlier_samples <- o.s
outlier_filter <- o.f# Add color sidebar
a <- 0.92
b <- 1
x <- tmp.sample1a[,6]
left.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
sidebar<-cbind(left.type, left.type, left.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
cbind(
cor.tmp[order(tmp.sample1a[,6]),
order(tmp.sample1a[,6])],
sidebar[order(tmp.sample1a[,6]),]),
col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Left Pos, z=(0.92,1)"), asp = 1)# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,7]
right.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
sidebar<-cbind(right.type, right.type, right.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
cbind(
cor.tmp[order(tmp.sample1a[,7]),
order(tmp.sample1a[,7])],
sidebar[order(tmp.sample1a[,7]),]),
col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Right Pos, z=(0.92,1)"), asp = 1)# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,8]
min.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
sidebar<-cbind(min.type, min.type, min.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
cbind(
cor.tmp[order(tmp.sample1a[,8]),
order(tmp.sample1a[,8])],
sidebar[order(tmp.sample1a[,8]),]),
col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Min Pos, z=(0.92,1)"), asp = 1)# Boxplot
plot(apply(cor.tmp[order(tmp.sample1a[,8]),order(tmp.sample1a[,8])],1,median))x <- apply(cor.tmp[order(tmp.sample1a[,8]),order(tmp.sample1a[,8])],1,median)
df <- data.frame(sample_id = names(x), cor_median = x, min_pos = sort(as.integer(tmp.sample1a[,8]))) %>%
left_join(y = pla_rppos_pass1a.0.order) %>%
mutate(min_pos = factor(min_pos))
df %>% ggplot(aes(x = min_pos, y = cor_median)) +
geom_boxplot(alpha = 0.5) +
geom_jitter(width = 0.2, height = 0, alpha = 0.8, color = 'forestgreen')# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,9]
sum.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
sidebar<-cbind(sum.type, sum.type, sum.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
cbind(
cor.tmp[order(tmp.sample1a[,9]),
order(tmp.sample1a[,9])],
sidebar[order(tmp.sample1a[,9]),]),
col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Block Length, z=(0.92,1)"), asp = 1)# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,10]
left.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
sidebar<-cbind(left.type, left.type, left.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
cbind(
cor.tmp[order(tmp.sample1a[,10]),
order(tmp.sample1a[,10])],
sidebar[order(tmp.sample1a[,10]),]),
col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Left Pos (Drift), z=(0.92,1)"), asp = 1)plot(apply(cor.tmp[order(tmp.sample1a[,10]),order(tmp.sample1a[,10])],1,median))plot(apply(cor.tmp[order(tmp.sample1a[,10]),order(tmp.sample1a[,10])],1,median), ylim = c(0.925,1))# Boxplot
x <- apply(cor.tmp[order(tmp.sample1a[,10]),order(tmp.sample1a[,10])],1,median)
df <- data.frame(sample_id = names(x), cor_median = x, left_pos = sort(as.integer(tmp.sample1a[,10]))) %>%
left_join(y = pla_rppos_pass1a.0.order) %>%
mutate(left_pos = factor(left_pos))
df %>% ggplot(aes(x = left_pos, y = cor_median)) +
geom_boxplot(alpha = 0.5) +
geom_jitter(width = 0.2, height = 0, alpha = 0.8, color = 'forestgreen')# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,11]
right.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
sidebar<-cbind(right.type, right.type, right.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
cbind(
cor.tmp[order(tmp.sample1a[,11]),
order(tmp.sample1a[,11])],
sidebar[order(tmp.sample1a[,11]),]),
col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Right Pos (Drift), z=(0.92,1)"), asp = 1)# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,12]
min.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
sidebar<-cbind(min.type, min.type, min.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
cbind(
cor.tmp[order(tmp.sample1a[,12]),
order(tmp.sample1a[,12])],
sidebar[order(tmp.sample1a[,12]),]),
col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Min Pos (Drift), z=(0.92,1)"), asp = 1)# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,13]
sum.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
sidebar<-cbind(sum.type, sum.type, sum.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
cbind(
cor.tmp[order(tmp.sample1a[,13]),
order(tmp.sample1a[,13])],
sidebar[order(tmp.sample1a[,13]),]),
col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Block Length (Drift), z=(0.92,1)"), asp = 1)plot(apply(cor.tmp[order(tmp.sample1a[,13]),order(tmp.sample1a[,13])],1,median))# Boxplot
x <- apply(cor.tmp[order(tmp.sample1a[,13]),order(tmp.sample1a[,13])],1,median)
df <- data.frame(sample_id = names(x), cor_median = x, drift_block_length = sort(as.integer(tmp.sample1a[,13]))) %>%
left_join(y = pla_rppos_pass1a.0.order) %>%
mutate(drift_block_length = factor(drift_block_length))
df %>% ggplot(aes(x = drift_block_length, y = cor_median)) +
geom_boxplot(alpha = 0.5) +
geom_jitter(width = 0.2, height = 0, alpha = 0.8, color = 'forestgreen')Note: Samples as columns Sidebar: p-values, t-values, and sig (binary) comparing normal test samples and outlier test samples
n.s <- row.names(pla1a.0.nr2)[!row.names(pla1a.0.nr2) %in% o.s]
# dim(pla1a.0.nr2[o.s,])
# dim(pla1a.0.nr2[n.s,])
hmo <- heatmap(pla1a.0.nr2)$colIndif(length(o.s) > 1){
i <- 1
df_t <- data.frame()
for(i in 1:ncol(pla1a.0.nr2)){
feat <- colnames(pla1a.0.nr2)[i]
tobj <- t.test(y = pla1a.0.nr2[n.s,i], x = pla1a.0.nr2[o.s,i])
df_t <- rbind(df_t, data.frame(feat, t_val = round(tobj$statistic,digits = 2),p_val = round(tobj$p.value, digits = 16)))
}
df_t <- df_t %>% arrange(t_val) %>%
mutate(sig = ifelse(p_val <= 0.05, 1, 0))
row.names(df_t) <- df_t$feat
# Add color sidebar for pval
a <- range(pla1a.0.nr2)[1]
b <- range(pla1a.0.nr2)[2]
x <- df_t[colnames(pla1a.0.nr2)[hmo], "p_val"]
p.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
x <- df_t[colnames(pla1a.0.nr2)[hmo], "t_val"]
t.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
x <- df_t[colnames(pla1a.0.nr2)[hmo], "sig"]
sig.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
sidebar<-rbind(p.type, p.type,
t.type, t.type,
sig.type, sig.type)
image(rbind(pla1a.0.nr2[order(tmp.sample1a[,2]),hmo],sidebar)
,col=redblue100,axes=F,main=glue("{TIS}2, f{P1}-HC, Run-Order"), asp = 1)
}else{
image(pla1a.0.nr2[order(tmp.sample1a[,2]),hmo],col=redblue100,axes=F,main=glue("{TIS}2, f{P1}-HC, Run-Order"), asp = 1)
}plot(apply(pla1a.0.nr2[order(tmp.sample1a[,2]),hmo],1,median))Note: pla1a.0.2 and pla1a.0.nr2 represent knn-imputed and log2 Note: Samples as columns
a <- range(pla1a.0.nr2)[1]
b <- range(pla1a.0.nr2)[2]
x <- tmp.sample1a[,3]
sex.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
x <- tmp.sample1a[,4]
control.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
x <- as.factor(tmp.sample1a[,5]) %>% as.numeric()
time.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
sidebar<-cbind(sex.type,sex.type,sex.type,sex.type,sex.type,
control.type,control.type, control.type, control.type, control.type,
time.type, time.type, time.type, time.type, time.type)
image(
cbind(pla1a.0.nr2[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]), heatmap(pla1a.0.nr2)$colInd ],
sidebar[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]),]),
col=redblue100,axes=F,main=glue("{TIS}2, f{P1}-HC, Sex-Control-Time"), asp =1)plot(apply(pla1a.0.nr2[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]),
heatmap(pla1a.0.nr2)$colInd],1,median))Strategy is not functionally programed (must revert log transformed back to linear gor pla1a.03a)
# pla1a.0.r3a<-pla1a.0.r3b<-pla1a.0.r3c<-pla1a.0.r3c2<-pla1a.0.r2b<-pla1a.0.r3d<-pla1a.0.r3d2<-pla1a.0.r2b2<-pla1a.0.r3d3 <-pla1a.0.2
pla1a.03a<-pla1a.03b<-pla1a.03c<-pla1a.03c2<-pla1a.02b<-pla1a.03d<-pla1a.03d2<-pla1a.02b2<-pla1a.03d3 <-pla1a.0.nr2
#pla1a.03c3<-pla1a.03b2<-pla1a.0.nr2tmp.s.median <- apply(pla1a.0.nr2,1, median)
tmp.s.mean <- apply(pla1a.0.nr2,1, mean)
plot(tmp.s.median,tmp.s.mean, asp = 1); abline(0,1)tmp.f.median <- apply(2^pla1a.0.nr2,2, median)
tmp.f.sd <- apply(2^pla1a.0.nr2,2, sd)
plot(y = tmp.f.sd, x = tmp.f.median,log="xy")tmp <- 2^pla1a.0.nr2[n.s,]
#dim(tmp)
tmp2.f.median <- apply(tmp,2, median)
tmp2.f.sd <- apply(tmp,2, sd)
plot(tmp.f.median,tmp2.f.median,log="xy")plot(tmp.f.sd,tmp2.f.sd,log="xy")for (i in 1:length(tmp.f.median)) {
pla1a.03a[,i]<-(2^pla1a.0.nr2[,i]- tmp.f.median[i])/ tmp.f.sd[i]
}
plot(2^pla1a.0.nr2[,1],pla1a.03a[,1]) #spot-check, verifiedNotice the drift downward
# Run Order; Original HC
###########################
hmo <- heatmap(pla1a.0.nr2)$colIndimage(pla1a.03a[order(tmp.sample1a[,2]), hmo],
col=redblue100,axes=F,main=glue("{TIS}3a, f{P1}-HC, Run-Order"), asp = 1)# Run Order; New HC
###########################
hmo <- heatmap(pla1a.03a)$colIndimage(pla1a.03a[order(tmp.sample1a[,2]), hmo],
col=redblue100,axes=F,main=glue("{TIS}3a, f{P1}-HC, Run-Order"), asp = 1)plot(apply(pla1a.03a[order(tmp.sample1a[,2]),
hmo],1,median)); abline(h = 0, col = 'blue')tmp.f.median <- apply(pla1a.0.nr2,2, median)
tmp.f.sd <- apply(pla1a.0.nr2,2, sd)
plot(tmp.f.median,tmp.f.sd, log = "xy")tmp <- pla1a.0.nr2[n.s,]
#dim(tmp)
tmp2.f.median <- apply(tmp,2, median)
tmp2.f.sd <- apply(tmp,2, sd)
plot(tmp.f.median,tmp2.f.median,log="xy")plot(tmp.f.sd,tmp2.f.sd,log="xy")for (i in 1:length(tmp.f.median)) {
pla1a.03b[,i]<-(pla1a.0.nr2[,i]- tmp.f.median[i])/ tmp.f.sd[i]
}
plot(pla1a.0.nr2[,1],pla1a.03b[,1]) #verified# Run Order; Original HC
###########################
hmo <- heatmap(pla1a.0.nr2)$colIndimage(pla1a.03b[order(tmp.sample1a[,2]), hmo],
col=redblue100,axes=F,main=glue("{TIS}3b, f{P1}-HC, Run-Order"), asp = 1)# Run Order; New HC
###########################
hmo <- heatmap(pla1a.03b)$colIndimage(pla1a.03b[order(tmp.sample1a[,2]), hmo],
col=redblue100,axes=F,main=glue("{TIS}3b, f{P1}-HC-NEW, Run-Order"), asp = 1)plot(apply(pla1a.03b[order(tmp.sample1a[,2]),
hmo],1,median))# Sex Control Time
a <- range(pla1a.03b)[1]
b <- range(pla1a.03b)[2]
x <- tmp.sample1a[,3]
sex.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
x <- tmp.sample1a[,4]
control.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
x <- as.factor(tmp.sample1a[,5]) %>% as.numeric()
time.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a)
sidebar<-cbind(sex.type,sex.type,sex.type,sex.type,sex.type,
control.type,control.type, control.type, control.type, control.type,
time.type, time.type, time.type, time.type, time.type)
image(
cbind(pla1a.03b[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]), heatmap(pla1a.0.nr2)$colInd ],
sidebar[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]),]),
col=redblue100,axes=F,main=glue("{TIS}3b, f{P1}-HC, Sex-Control-Time"), asp = 1)plot(pla1a.03b[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]),])Outlier sample has increased sd
boxplot(data.frame(t(pla1a.03b)))plot(apply(pla1a.0.nr2,1,median), apply(pla1a.03b,1,median))# remove outliers and determine if there is significant variablility in sample medians
glue("Outlier samples removed:")## Outlier samples removed:
boxplot(data.frame(t(pla1a.03b[n.s,])))# # log 2
# pla_rppos_pass1a_nr1 <- pla1a.0.nr1
# pla_rppos_pass1a_1 <- pla1a.0.1
#
# # Test only
# #############
# # log2, imputed (test only)
# pla_rppos_pass1a_nr2 <- pla1a.0.nr2
# # linear, median-centered, sd-scaled
# pla_rppos_pass1a_3a <- pla1a.03a
# # log2, median-centered, sd-scaled
# pla_rppos_pass1a_3b <- pla1a.03b
feature_impute = names(feature_impute)
# The processing decisions
pla1a.0_vars <- data.frame(ds, site, tech, tis, NR, N, P, neg_vals, zero_vals, feature_na_filter,
P1, NR1, N1, paste(feature_impute, collapse = '; '), na_vals_impute, knn_k,
run_var, outlier_sample_n = length(o.s), outlier_samples = paste(o.s, collapse = '; '), outlier_filter = o.f,
internal_standards = paste(is, collapse = '; '),
comments = "No obvious sample outliers, data looks clean")
# visualize the processing decisions
if(knit_time){
pla1a.0_vars %>%
knitr::kable(format = "html") %>%
scroll_box(width = "100%", height = "100%")
}| ds | site | tech | tis | NR | N | P | neg_vals | zero_vals | feature_na_filter | P1 | NR1 | N1 | paste.feature_impute..collapse…….. | na_vals_impute | knn_k | run_var | outlier_sample_n | outlier_samples | outlier_filter | internal_standards | comments |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| pass1a | umich | rppos | pla | 101 | 77 | 170 | 0 | 0 | 10 | 170 | 101 | 77 | 0 | 10 | 0 | 1 | 90122013104 | 0.96 | valine-13C5 [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; lysine-13C6 [iSTD]; methionine-13C5 [iSTD]; Car-D9 [iSTD]; phenylalanine-13C9 [iSTD]; arginine-13C6 [iSTD]; glucose-13C6 [iSTD]; tyrosine-13C9 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Car(2:0)-D3 [iSTD]; Car(3:0)-D3 [iSTD]; Car(4:0)-D3 [iSTD]; Car(8:0)-D3 [iSTD]; fructose 1,6-bisphosphate-13C6 [iSTD]; Gibberellic acid [iSTD]; Car(16:0)-D3 [iSTD]; 24-epi-Brassinolide [iSTD] | No obvious sample outliers, data looks clean |
# save the data
save(pla1a.0.nr1, pla1a.0.1, pla1a.0.nr2, pla1a.03a, pla1a.03b, pla1a.0_vars,
file = glue("{WD}/data/UM-rppos/UM_rppos_processed_pla1a.0.0.RData"))warnings()
session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.0.2 (2020-06-22)
## os macOS 10.16
## system x86_64, darwin17.0
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz America/Toronto
## date 2021-12-11
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date lib source
## assertthat 0.2.1 2019-03-21 [2] CRAN (R 4.0.2)
## backports 1.2.1 2020-12-09 [2] CRAN (R 4.0.2)
## broom 0.7.8 2021-06-24 [1] CRAN (R 4.0.2)
## cachem 1.0.3 2021-02-04 [2] CRAN (R 4.0.2)
## callr 3.7.0 2021-04-20 [1] CRAN (R 4.0.2)
## cellranger 1.1.0 2016-07-27 [2] CRAN (R 4.0.2)
## cli 3.0.1 2021-07-17 [1] CRAN (R 4.0.2)
## coda 0.19-4 2020-09-30 [1] CRAN (R 4.0.2)
## codetools 0.2-18 2020-11-04 [2] CRAN (R 4.0.2)
## colorspace 2.0-2 2021-06-24 [1] CRAN (R 4.0.2)
## crayon 1.4.1 2021-02-08 [2] CRAN (R 4.0.2)
## curl 4.3.2 2021-06-23 [1] CRAN (R 4.0.2)
## data.table 1.14.0 2021-02-21 [1] CRAN (R 4.0.2)
## DBI 1.1.1 2021-01-15 [2] CRAN (R 4.0.2)
## dbplyr 2.1.1 2021-04-06 [1] CRAN (R 4.0.2)
## desc 1.3.0 2021-03-05 [1] CRAN (R 4.0.2)
## devtools * 2.4.2 2021-06-07 [1] CRAN (R 4.0.2)
## digest 0.6.27 2020-10-24 [2] CRAN (R 4.0.2)
## dplyr * 1.0.7 2021-06-18 [1] CRAN (R 4.0.2)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.0.2)
## evaluate 0.14 2019-05-28 [2] CRAN (R 4.0.1)
## fansi 0.5.0 2021-05-25 [1] CRAN (R 4.0.2)
## farver 2.1.0 2021-02-28 [1] CRAN (R 4.0.2)
## fastmap 1.1.0 2021-01-25 [2] CRAN (R 4.0.2)
## forcats * 0.5.1 2021-01-27 [2] CRAN (R 4.0.2)
## fs 1.5.0 2020-07-31 [2] CRAN (R 4.0.2)
## generics 0.1.0 2020-10-31 [2] CRAN (R 4.0.2)
## ggfittext 0.9.1 2021-01-30 [2] CRAN (R 4.0.2)
## ggplot2 * 3.3.5 2021-06-25 [1] CRAN (R 4.0.2)
## glue * 1.4.2 2020-08-27 [2] CRAN (R 4.0.2)
## gridExtra 2.3 2017-09-09 [2] CRAN (R 4.0.2)
## gtable 0.3.0 2019-03-25 [2] CRAN (R 4.0.2)
## haven 2.3.1 2020-06-01 [2] CRAN (R 4.0.2)
## highr 0.9 2021-04-16 [1] CRAN (R 4.0.2)
## hms 1.1.0 2021-05-17 [1] CRAN (R 4.0.2)
## htmltools 0.5.1.1 2021-01-22 [2] CRAN (R 4.0.2)
## httr 1.4.2 2020-07-20 [2] CRAN (R 4.0.2)
## impute * 1.62.0 2020-04-27 [1] Bioconductor
## inline 0.3.19 2021-05-31 [1] CRAN (R 4.0.2)
## inspectdf 0.0.11 2021-04-02 [1] CRAN (R 4.0.2)
## jsonlite 1.7.2 2020-12-09 [2] CRAN (R 4.0.2)
## kableExtra * 1.3.1 2020-10-22 [2] CRAN (R 4.0.2)
## knitr 1.33 2021-04-24 [1] CRAN (R 4.0.2)
## labeling 0.4.2 2020-10-20 [2] CRAN (R 4.0.2)
## lattice 0.20-41 2020-04-02 [2] CRAN (R 4.0.2)
## lifecycle 1.0.0 2021-02-15 [1] CRAN (R 4.0.2)
## loo 2.4.1 2020-12-09 [1] CRAN (R 4.0.2)
## lubridate 1.7.10 2021-02-26 [1] CRAN (R 4.0.2)
## magrittr 2.0.1 2020-11-17 [2] CRAN (R 4.0.2)
## MASS 7.3-53 2020-09-09 [2] CRAN (R 4.0.2)
## matrixStats 0.60.0 2021-07-26 [1] CRAN (R 4.0.2)
## memoise 2.0.0 2021-01-26 [2] CRAN (R 4.0.2)
## modelr 0.1.8 2020-05-19 [2] CRAN (R 4.0.2)
## MotrpacBicQC * 0.5.2 2021-06-25 [1] Github (MoTrPAC/MotrpacBicQC@43fb293)
## munsell 0.5.0 2018-06-12 [2] CRAN (R 4.0.2)
## mvtnorm 1.1-2 2021-06-07 [1] CRAN (R 4.0.2)
## naniar 0.6.1 2021-05-14 [1] CRAN (R 4.0.2)
## pillar 1.6.2 2021-07-29 [1] CRAN (R 4.0.2)
## pkgbuild 1.2.0 2020-12-15 [2] CRAN (R 4.0.2)
## pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 4.0.2)
## pkgload 1.2.1 2021-04-06 [1] CRAN (R 4.0.2)
## plyr 1.8.6 2020-03-03 [2] CRAN (R 4.0.2)
## prettyunits 1.1.1 2020-01-24 [2] CRAN (R 4.0.2)
## processx 3.5.2 2021-04-30 [1] CRAN (R 4.0.2)
## progress 1.2.2 2019-05-16 [2] CRAN (R 4.0.2)
## ps 1.6.0 2021-02-28 [1] CRAN (R 4.0.2)
## purrr * 0.3.4 2020-04-17 [2] CRAN (R 4.0.2)
## R6 2.5.0 2020-10-28 [2] CRAN (R 4.0.2)
## Rcpp 1.0.7 2021-07-07 [1] CRAN (R 4.0.2)
## RcppParallel 5.1.4 2021-05-04 [1] CRAN (R 4.0.2)
## readr * 1.4.0 2020-10-05 [2] CRAN (R 4.0.2)
## readxl 1.3.1 2019-03-13 [2] CRAN (R 4.0.2)
## remotes 2.4.0 2021-06-02 [1] CRAN (R 4.0.2)
## reprex 2.0.0 2021-04-02 [1] CRAN (R 4.0.2)
## reshape2 1.4.4 2020-04-09 [2] CRAN (R 4.0.2)
## rethinking * 2.13 2021-08-13 [1] Github (rmcelreath/rethinking@2acf2fd)
## rlang 0.4.11 2021-04-30 [1] CRAN (R 4.0.2)
## rmarkdown 2.6 2020-12-14 [2] CRAN (R 4.0.2)
## rprojroot 2.0.2 2020-11-15 [2] CRAN (R 4.0.2)
## rstan * 2.21.2 2020-07-27 [1] CRAN (R 4.0.2)
## rstudioapi 0.13 2020-11-12 [2] CRAN (R 4.0.2)
## rvest 1.0.0 2021-03-09 [1] CRAN (R 4.0.2)
## scales 1.1.1 2020-05-11 [2] CRAN (R 4.0.2)
## sessioninfo 1.1.1 2018-11-05 [2] CRAN (R 4.0.2)
## shape 1.4.6 2021-05-19 [1] CRAN (R 4.0.2)
## StanHeaders * 2.21.0-7 2020-12-17 [1] CRAN (R 4.0.2)
## stringi 1.7.3 2021-07-16 [1] CRAN (R 4.0.2)
## stringr * 1.4.0 2019-02-10 [2] CRAN (R 4.0.2)
## testthat 3.0.4 2021-07-01 [1] CRAN (R 4.0.2)
## tibble * 3.1.3 2021-07-23 [1] CRAN (R 4.0.2)
## tidyr * 1.1.3 2021-03-03 [1] CRAN (R 4.0.2)
## tidyselect 1.1.1 2021-04-30 [1] CRAN (R 4.0.2)
## tidyverse * 1.3.1 2021-04-15 [1] CRAN (R 4.0.2)
## usethis * 2.0.1 2021-02-10 [1] CRAN (R 4.0.2)
## utf8 1.2.2 2021-07-24 [1] CRAN (R 4.0.2)
## V8 3.4.2 2021-05-01 [1] CRAN (R 4.0.2)
## vctrs 0.3.8 2021-04-29 [1] CRAN (R 4.0.2)
## viridisLite 0.4.0 2021-04-13 [1] CRAN (R 4.0.2)
## visdat 0.5.3 2019-02-15 [2] CRAN (R 4.0.2)
## webshot 0.5.2 2019-11-22 [2] CRAN (R 4.0.2)
## withr 2.4.2 2021-04-18 [1] CRAN (R 4.0.2)
## xfun 0.24 2021-06-15 [1] CRAN (R 4.0.2)
## xml2 1.3.2 2020-04-23 [2] CRAN (R 4.0.2)
## yaml 2.2.1 2020-02-01 [2] CRAN (R 4.0.2)
##
## [1] /Users/Alec/Library/R/4.0/library
## [2] /Library/Frameworks/R.framework/Versions/4.0/Resources/library